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The Ising lattice gas, with its well known equilibrium properties, displays a number of surprising 
phenomena when driven into non-equilibrium steady states. We study such a model with anisotropic 
interparticle interactions ( J« J± ), using both Monte Carlo simulations and high temperature 
series techniques. Under saturation drive, the shift in the transition temperature can be both 
positive and negative, depending on the ratio Jy/Jx! For finite drives, both first and second order 
transitions are observed. Some aspects of the phase diagram can be predicted by investigating the 
two point correlation function at the first non-trivial order of a high temperature series expansion. 



I. INTRODUCTION 



Our understanding of collective behavior in many-body systems in, or near, thermal equilibrium has improved 
significantly in the past century. By contrast, only recently has there been much attention focused on systems far- 
from-equilibrium, even for those in steady states. Abundant in nature, such systems exist in an enormous variety 
of states, most of which cannot be predicted using the framework of equilibrium statistical mechanics. At present, 
there is no simple formulation which is generally applicable. One way to make progress into this vast realm of non- 
equilibrium physics is to study simple models, especially those with well understood equilibrium properties. This is 
the chief motivation of Katz, et.al., who introduced a seemingly trivial modification to the Ising lattice gas so 
that it is driven into non- equilibrium steady states. The other motivation is phase transitions in fast ionic conductors 
in an external field, though their model predictions have not been compared directly with these physical systems so 
far. 

This model consists of a <i-dimcnsional hyper-cubic lattice with each site being empty or occupied by a single 
particle. Its dynamics is particle hopping to nearest neighbor empty sites, with a rate controlled by the energetics of 
the Ising Hamiltonian Tt, as well as a bias in one direction so as to describe the effect of a uniform, DC "electric" field 
E. To connect with the equilibrium system, rates are chosen to satisfy "local" detailed balance consistent with being 
in contact with a thermal bath at temperature T . In particular, by setting E = 0, we simply retrieve the canonical 
exp (— H/huT) as the stationary distribution. Since its inception, many unexpected properties have been discovered, 
for both the original model and its many variants |3| . Though some of the remarkable behaviors are reasonably well 
understood, many remain unexplained. An example is the basic question: why should T C (E), the critical temperature 
for the attractive (ferromagnetic) case, increase with El Simulations (using Metropolis rates for d — 2) show that 
T c (oo) is about 40% above the Onsager temperature 0! Indeed, one can easily argue in favor of a lowering of T c , as 
follows. Since large fields would dominate over the inter-particle attraction (for hops along E), the system is effectively 
subjected to an extra noise, so that a lower T is needed for clustering and order. 

With further explorations of this system, simple arguments in favor of an increased T c began to emerge [g| . However, 
there is still no compelling reason, so far, to accept or reject the contradictory arguments. In other words, we still 
have no intuitive picture which can guide us to reliable expectations. This study is a continuation of the exploration 
of similar systems in an effort to find the underlying mechanisms which give rise to the novel phase diagrams. 

Specifically, we consider Ising models with anisotropic interactions, i.e., inter-particle attraction along and transverse 
to the drive being unequal (J|| ^ J±)- Such a model was previously studied, but only in the fast rate limit [||. 
Unfortunately, though very interesting in its own right, this limit is too singular to provide us with much insight into 
the original lattice gas. For example, no trace of E remains, yet a phase transition prevails even for J± = 0! In the 
following section, we describe our model, restricting our attention to d = 2 only. Simulations lead to more surprises, 
the details of which will be presented in Section 3. The next section is devoted to a theoretical study, using high 
temperature series expansion techniques on the two-point correlation function J9|,[l0|]. A final section is devoted to a 
summary, conclusions and outlook. 
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II. DRIVEN LATTICE GAS WITH ANISOTROPIC INTERACTIONS 



Our system consists of fully periodic L x M square lattices, with the sites labeled by x = (x, y), where x and y are 
integers modulo L and M , respectively. Each site may be empty or occupied by a particle, so that a configuration of 
the system is specified by the set of occupation numbers {n(x)}, where n = or 1. This model is simply related to the 
original Ising one for spins through the relation s = 2n — 1, which takes on values ±1. Since our interest is a system 
of particles, ^ x n(x) is fixed. For simplicity, we study only half-filled systems, i.e., ^ n = LM/2 (or s — 0). Next, 
we endow the particles with anisotropic nearest neighbor attraction and write the Hamiltonian as 

H = -4J X ^2 n ( x > y) n ( x + 1) y) - 4J y ^2 n ( x ' v) n ( x > v + l ) (!) 

with J x ,J y > 0. The factor of 4 is included so that the Hamiltonian is equivalent (given the constraint Yl s = 0) 
to the standard Ising form, —J^ss'. For convenience, we introduce anisotropy via a "dimensionless" parameter a, 
through 

Jx = J/ca and J y = J a (2) 

so that J carries the "overall" strength of the interactions. When this system reaches equilibrium with a thermal 
bath at temperature T, it displays well known properties in the thermodynamic limit 0-03). The most prominent 
behavior is a second order transition, from a disordered homogeneous phase to a phase segregated state, when T is 
lowered through the Onsager temperature T c (a). Expressed in units of J/fcs, where ks is Boltzmann's constant, 
T c (a) can be obtained from the equation 

(l+e Q ) = 2 (3) 

for the quantity e = exp (— 2 J/ksT c (a)). Specifically, T c (l) ~ 2.269J/fcs, which serves as the unit for all the 
temperatures quoted in this paper. Using this unit, a plot of T c (a) is provided in Fig. 1. Note that, due to the 
conservation law, the ordered state will be a strip spanning the system aligned with the x- or ?/-axis. These states 
will be denotes by H (horizontal) and V (vertical), respectively. Since the equilibrium state will be the one with the 
lowest interfacial free energy, the aspect ratio of the system, L/M, will play a role as well. For L/M = 1, clearly 
a = 1 marks the "phase boundary" between the two states, also shown in Fig. 2. We expect, of course, that the H-V 
boundary is associated with a first order transition. For general L/M, this boundary is located at that a(T) which 
satisfies a y (a,T)/a x (a,T) = L/M, where a x is the surface tension for the horizontal interface, etc. 

In computer studies, the simulation of the lattice gas, in equilibrium with a bath at temperature T, would rely on 
"spin-exchange" dynamics |l5[ ], so as to respect particle conservation. In other words, particles are allowed to hop to 
nearest neighbor holes with a rate obeying detailed balance appropriate for T. A favorite rate is due to Metropolis 
|l6|| , where jump rates are given by min{l, exp (— ATl/ksT)}, where ATI is the change in energy due to the hop. In 
the generalization to a driven lattice gas, we follow Katz, et. al. and incorporate the "electric" field (aligned with the 
y-axis) in the standard way, i.e., by adding ±E to AH for hops against /along the field Jl|,||. Our goal is to map out 
the phase diagram in the T-a-E space. In the following, we will report the main results, denoting the line of second 
(first) order transitions by T c (a, E) (ai(T, E)) and quoting E in units of J. Details will be published elsewhere |p7| . 



III. SIMULATION METHODS AND RESULTS 



Most of our runs involve lattices with L — M — 30. For definiteness, we aligned our drive so that particles are 
biased to hop "downwards" (— y direction). Starting with one of three initial configurations of a half-filled lattice, 
the system is evolved in the standard way. In a Monte Carlo step (MCS), 2LM bonds (nearest neighbor pairs) 
are chosen at random. If the bond consists of a particle hole pair, then the pair is exchanged with probability 
min{l, exp (— [ATI + ESy] /fc^T)}, where Sy is the change in the y-coordinate (mod M) of the particle. Most runs are 
400K MCS long. The initial configurations are either random or fully ordered (in H or V). Typically, we discard the 
first 100K MCS, allowing the system to settle into the steady state. Then, every 200 MCS, we measure the Fourier 
transform of the particle density n(x) 



i( k ,P) = -^=^2n(x,y) exp 



. kx py 



(4) 



where k,p are integers. Averaging over the rest of the run (denoted by ()), we obtain the structure factors 
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S(k,p) = (|n(fc,p)| 2 ). 



(5) 



In the disordered state, they are of order unity and convey information about two-particle correlations. Below the 
transition, they are sensitive to the H- or V-ordered state, in the sense that, e.g., 5(0, 1) or 5(1,0) will be O(LM), 
and carry information about the densities of phase segregated states. For example, for a completely ordered V state, 
we obtain 5(0, 1) = while 5(1,0) = (2M/L) [1 — cos(27r/L)]~ 1 , which approaches LAl/ir 2 in the thermodynamic 
limit. Of course, in an ordered H state, these two 5's are reversed. In this sense, 5(0, 1) and 5(1,0) may be viewed 
as order parameters In addition to the mean of \n\ 2 , we measure the variance, (A5) = (|n| 4 ) — 5 2 , which plays a 
role similar to that of the susceptibility in equilibrium systems. For those 5's which are 0(1), the associated A5's are 
expected to be just 5 itself [Q. Deep in the ordered phases, though some 5's will be O(LM), the A5's are expected 
to be still 0(1). However, near a second order transition, the A5's should diverge in the thermodynamic limit. We 
have also exploited other initial configurations, different sizes and rectangular systems. These will be discussed in the 
following sub-sections. 

A. Second order transitions with saturation drive 

In the first part of our study, we restrict ourselves to saturation drive ( E = 50 here) and to only one system 
size/geometry: 30 x 30. Since the field is large enough to overcome all levels of particle attraction, it is hardly 
surprising that the ordered state can only be V. Measuring the structure factors, we confirm this expectation, in that 
only S(k, 0) (k odd) grow substantially as T is lowered. Specifically, starting with random configurations, we perform 
400 K MCS runs for a = |, i, |, 1, |, 2, 3 and 0.5 < T < 3.0. We used a range of step sizes in T, the smallest being 
0.025, especially when large variations in measured quantities were encountered. Identifying the transition through the 
peak of (A5(l,0)) , we find the line of second order transitions T c (a,oo) to be a monotonically decreasing function 
of a, as shown in Fig. 1. This behavior is contrary to our previous conjecture j^], showing that our intuitive picture, 
based on the competition of long and short ranged correlations, is still far from being reliable. More surprisingly, the 
critical temperature of the driven system is not always greater than that for the equilibrium system! Note that, in 
Fig. 1., T c (a, 0) is the theoretical result for an infinite lattice, while T c (a,oo) is found in the 30 x 30 system. (The 
effects of finite size will be discussed in a later sub-section.) A further distinction is that, while the driven system 
orders only into the V state, the equilibrium system may order into H or V depending on a < 1 or > 1. 

We should caution that the uncertainties associated with T c (a, oo) are much larger than the step size of 0.025. In 
particular, for T greater than the quoted T c (a,<x>), 5(1,0) does not decrease in a simple way, especially for systems 
with large a. In these cases, we frequently observe breakup of a single strip into two-strip states. As a result, as T is 
raised through this T c (a, oo) , not only does (A5(l, 0)) 2 becomes large, 5(2, 0) also increases by as much as an order 
of magnitude. In this sense, it is possible to regard the quoted values as a lower bound for the critical temperature. 
Typically, at temperatures ~ 30% higher than T c (a, oo), we can be quite certain that the system is disordered. To 
be more confident and precise about critical temperatures, we must use more sophisticated methods, e.g., histograms 
of various 5's, finite size scaling, etc. Q 

B. First and second order transitions with small drives 

Next, we turn to systems with small drives, searching for the remnants of the H phase. Keeping the same system 
size (30 x 30), we chose E = 1,2, and 5, so that there is a range of values for the ratio E/J± = olE. Expecting 
continuous transitions in the neighborhood of the equilibrium line, we performed simulations in the manner described 
above. As displayed in Fig. 2, the behavior of T c (a,E) is quite complex. For a — 3, T c first increases and then 
decreases, so as to match with the low value of T c (3,oo) found above. For a = 1, the known monotonic increase of 
T c is confirmed. In its neig hborhood, both T C (±,E) and T C (\,E) appear to be monotonic in E. In the latter case, 
the ordered phase remains V, despite the fact a < 1. This reflects the known effect of the drive: setting up positive 
long range correlations along the drive and so, favoring strips aligned vertically. For a < 3/4 and the E's we have 
chosen, the systems orders into V or H depending on the drive. For those cases where the ordered phase is H, it is 
completely unexpected that T c (a,E) is still typically greater than the equilibrium T c (a,0)\ Since the drive is known 
to reduce correlations in the transverse direction, we must find a new argument for how ordering (into H) could occur 
at a higher temperature. Remarkably, T c (a,E) suffers a dip at the bi-critical point. As a result, for a fixed a, it is 
seen first to increase with E (while ordering into H), then to decrease to the minimum value T c (a, Eb) when the drive 
reaches Eg, a strength just strong enough to change the ordering to V, and finally, to increase again with stronger 
drives. Given our accuracy, T c (cx,Eb) seems to be the same as T c (a,0), though we have no reason to believe that 
they are indeed equal! 
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Figure 1. Critical temperatures vs. £na for E = 5(D) and 50(»). Both ordered states are V. 
The exact values for E — are labelled by o's; the dashed line is the boundary between H- and V-states. 
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Figure 2. Phase Diagram for E — (solid and dashed lines), 1 (+), and 2 (k). 

Points associated with first order transitions are joined by a dotted line. 
The disordered region is labelled by D; the ordered ones by H and V (See text). 
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Turning to the line of discontinuous transitions, we modified our methods. Since it is expected to be more or less 
aligned with the T axis, we sweep in a with fixed T and E. The critical values are denoted by a\(T,E), shown as 
points joined by dotted lines in Fig. 2. 

For E = 1,2 and T as low as 0.9, we are able to look for hysteresis in such sweeps. In both cases, a is increased 
and decreased by steps of 0.025 and the system evolved for 200K MCS at each step. The ratio 

5(0,1) -5(1,0) 
5(0,1) + 5(1,0) 

is monitored and is seen to jump from 1 to -1 and vice versa. a\(T,E) is found by averaging the values of a where 
the jump occurs. 

For lower T, the mctastable life times are too long and hysteresis loops become too large to be reliable. Thus, we 
resort to another method, which, to our knowledge, is entirely new. Though it does not necessarily identify the loci 
of first order transitions, we argue that it should provide a good indication. This approach relies on a conjecture of 
the saddle point configuration, namely, an "X" (Fig. 3a), which is chosen via symmetry considerations. Given the 
anisotropies in the system (due to both a and E), we concede that the saddle point may be a configuration with less 
symmetry, e.g., an "L" (Fig. 3b). Nevertheless, we believe that the "X" should be a good starting point. First, we 
carry out 400K MCS runs with V and/or H as the initial configuration in this region of phase space. We chose only 
three values of T : 0.65, 0.70,and 0.75. For E = 1, we focused on the values a = 0.65, 0.70, 0.75, and 0.80. For E = 2, 
we used a = 0.40, 0.45, 0.50, and 0.55. Since the system remains in the initial state, sharp distributions of the structure 
factors (5(1, 0), 5(2, 0), 5(0, 1), and 5(0, 2)) emerge, providing us with good averages and standard deviations. Next, 
we carried out 200 independent runs for each of the above (T,a,E), starting with the "X" configuration. The above 
5's are measured every 200 MCS. The runs are terminated when 10 consecutive measurements (i.e., over a period of 
2 K MCS) of each 5 fall within 3 standard deviations of the averages. The fraction of runs which terminate in the V 
state is recorded and found to increase monotonically with a (with fixed T, E). The points where this fraction reaches 
1/2 are considered (part of) the line of first order transitions from H to V. (See Fig. 2.) 




Finally, for E = 5, all ordered states we have observed are V. Thus, we never explored the line of first order 
transitions. To keep Fig. 2 relatively clear, we display this line of continuous transitions in Fig.l instead. It is possible 
that first order transitions occur just outside the range of a's we used. Indeed, the slight dip of T c (a, 5) at a = 1/3 is 
actually due to behavior which is more complicated than that in systems with higher a. For T > 1.6, 5(1, 0) < 5(0, 1), 
even though the system actually orders into V for T < 1.5, i.e., 5(1,0) S> 5(0,1). In between, the fluctuations of 
both structure factors are comparably large. In this sense, it is possible that (a = 1/3, T ~ 1.55) is very close to 
the bi-critical point. To explore further, we must set a < 1/3. However, we believe that, with J±/J\\ > 10, reliable 
conclusions cannot be drawn from simulations with a lattice as small as 30 x 30 . 
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C. Effects of varying system sizes and aspect ratios 



The above data were collected on a single system size: 30 x 30. Of course, true phase transitions accompanied 
by thermodynamic singularities are properties of infinite systems. Any conclusions would rely on either simple 
extrapolations or more sophisticated methods of finite size scaling. With our limitations of computation power, we 
are able to make only cursory explorations of the effects of finite size. In particular, we studied 60 x 60 and 90 x 90 
systems, only for E = 50 and the extreme a's. The transition temperatures are summarized in Table 1. Similar to 
the equilibrium cases, the transition temperatures tend to rise with larger sizes. Although there are some significant 
increases, we can safely conclude, from the trends displayed here, that the general features of Figs.l and 2 will survive 
the thermodynamic limit. 

Table 1. Transition temperatures for various L and a. 



L \ a 


1/3 


1/2 


2 


3 


30 


2.25 


1.75 


0.97 


0.67 


60 


2.40 


1.85 


1.05 


0.77 


90 


2.45 


1.90 


1.05 


0.80 



Finally, we also investigated the effects of having rectangular systems. By varying the aspect ratio, we can insure 
that the ground state of the equilibrium system is in the V- or the H-state, with the crossover point at M/L = a 2 . 
We could study systems with these aspect ratios. On the other hand, at any finite temperature, the "crossover" M/L 
is given by the ratio of surface energies associated with the interfaces |19[ aligned along the two axes, i.e., 



M _ 2Ka + £n(temh(K/a)) 
T ~ 2K/a + £n (tanh (Ka)) 



K = J/k B T 



(6) 



In other words, for a given aspect ratio, this equation provides the phase boundary between H- and V-states. By 
duality, this ratio is also the one for the correlation lengths associated with two-point correlations along the two axes 
pot in the disordered phase. Thus, another possibility is to study systems which, near criticality, behave isotropically 
when the rectangular systems are rescaled to be square ones. Since our interest is mainly in disorder-order transitions, 
we adopted the latter choice here. Evaluating (Q) near T c , we arrived at lattice sizes 25 x 36, 20 x 48, and 15 x 64 
(for a — 4/3, 2, and 3 respectively), in order to compare with the 30 x 30 case at a — 1. By symmetry, systems with 
the opposite aspect ratio are used for a < 1. Keeping in mind the previous discussion concerning the uncertainties 
associated with T c , we display the results in Fig. 4. Note that the equilibrium line is not shown, since it is so close 
to the E = 1 points that confusion may arise. It is unclear if this insensitivity of T c (a, 1) to the drive is significant 
or not. Note that for a = 1/3 and 1/2, the system ordered into an H-state. For E = 2, the lowering of T c (l/3, 2) 
is curious. This point also coincides with ordering into an H-state. Comparing with Fig. 2, we see that there are 
two main effects due to rectangular lattices: (i) displacement of the first order line to smaller a's and (ii) lowering 
of the transition temperatures into the H-state. The lack of a dip in the E = 5 case is consistent with the bi-critical 
point being at a < 1/3. Apart from these features, system geometry seems to play no discernible role, especially 
for large drives. That the system subjected to small E orders into H-states may be related to the fact that, at low 
temperatures, the H-state is favored in the equilibrium systems according to eqn.(|^). Clearly, there remains a large 
gap between our understanding of these systems and the rich phenomena displayed. 



IV. TWO POINT CORRELATION IN A HIGH TEMPERATURE EXPANSION 



Though the driven lattice gas model was introduced 15 years ago, there is still no reliable means to predict the 
general features of the phase diagram. The dynamic mean-field approach |2l| , while quantitatively more satisfying, 
is so labor-intensive that it provides us with little insight on why critical temperatures shifts to higher values. This 
method can certainly be extended to our anisotropic system. Another route to an estimate of the shift of T c is 
based on a recent study of the two-point correlation function G(x) in a "high temperature" expansion []l0|] , which is 
actually an expansion in small J or K = J/ksT. In this approach, an approximate equation for G(x), first derived 
in H is solved exactly. The resultant Fourier transform is just the theoretical structure factor: 5'(k) which not only 
displays the well understood discontinuity singularity at k = [Q, but also may be exploited to estimate T C (E) ]To| |. 
The results are less accurate than those from dynamic mean-field theory, since the latter is the generalization of the 
Bethe-Peierls approximation for equilibrium Ising models. However, its implementation is much simpler. Here, we 
provide a brief presentation of the generalization to the anisotropic model. 
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Figure 4. Transition temperatures for rectangular systems as a function of £na. 
Different drives are shown as E = 1(+), 2(a), 5(D), and 50(»). 
Transitions into the V-state are joined by lines. 



With reference to pJTc|] , where the derivation and the nature of the approximation can be found, we simply quote 
the equation for G : 

d t G(0, 1) = 2[G(1, 1) + G(-l, 1) - 2G(0, 1)] + (1 + e)[G(0, 2) - G(0, 1)] + 2K(2 + e)a 

d t G(l, 0) = 2[G(2, 0) - G(l, 0)] + (1 + e)[G(l, 1) + G(l, -1) - 2G(1, 0)] + 2K(1 + 2e)/a 

d t G(l, 1) = 2[G(2, 1) + G(0, 1) - 2G(1, 1)] + (1 + e)[G(l, 2) + G(l, 0) - 2G(1, 1)] - 2K{a + e/a) (7) 

d t G(0, 2) = 2[G(1, 2) + G(-l, 2) - 2G(0, 2)] + (1 + e)[G(0, 3) + G(0, 1) - 2G(0, 2)] - 2eKa 

d t G(2, 0) = 2[G(3, 0) + G(l, 0) - 2G(2, 0)] + (1 + e)[G(2, 1) + G(2, -1) - 2G(2, 0)] - 2K/a 

and, for all other non-zero x, y: 

d t G(x, y) = (l + e)[G(x, y + l) + G(x, y - 1) - 2G(x, y)] + 2{G{x + l,y) + G(x - l,y) - 2G(x, y)] (8) 

Here, 

e = exp(-E/k B T) 

with the understanding that these equations are valid to first order in K. Setting the left hand side of (Q) to zero, 
the linear equations for the steady state G*(x,y) and its transform S*(k,p) can be solved, with the proviso that, at 
the zeroth order in K, 

G*(0,0) = 1 

G*(x,y) = x,y^0, 

and 

S*(k,p) = 1. 
7 



Since a does not affect the structure of these equations and enters only through the inhomogeneous terms, we can 
follow the analysis in |]lof closely and obtain the equation for T c (a, E). The result can be succinctly summarized by 

T c (a,E) _ 1 /(e) 1 

T C (1,E) l + /(e) a+ l + f(e)a W 

The structure of this ratio is necessarily of this form, since T c (a,E) is linear in (Ju , J±)=J(a, 1 / a) and, of course, 
the ratio is unity when a = 1. Also, since this approach reduces to the simple mean-field result [[24[ , a + 1/a, for the 
equilibrium case (E = 0), we know /(l) = 1. Unfortunately, the expression for / 

j,,, _ -gii + 2(1 - e) [R10R12 - 

i?n + 2(1 - e) [i?i i?2i — -R20-R11] 

where 

R^'TdMp (l-cos^l-cospT- 
J (27r) 2 i_ T 2(l + e)(l-cosA:)+4(l-cosp) 

is not transparent enough to provide much insight into shifts in T c . For the entire range < e < 1, we find that 
1 < / < 2, so that both coefficients in (||) are positive. As a result, T c (a) is not a monotonically decreasing function of 
a, as found in simulations for large E. Perhaps we should not be too surprised, since T c (a, E)/T c (l, E) can be only a 
linear combination of a and 1/a. Nevertheless, some features are qualitatively reproduced, such as T c (ao < 1) being 
greater than T c (l/ao). Finally, noticing that these results are based on an infinite lattice, we extended them to finite 
lattices. This generalization consists mainly of replacing the integrals in ( |ll| ) by finite sums over, e.g., k = 2nk/L. 
No significant changes are found down to L = 30, so that we did not extend this investigation to L 7^ M . Again, it is 
remarkable that qualitative features, such as the increase of T c with L, are reproduced. 

A more serious disadvantage of this approach is that it gives no information for E = O(J). If the "high temperature" 
expansion is also analytic in E, then the only way E can appear in T c (a,E) is through 0(E 2 ). For reasons beyond 
the scope of this paper, a second order computation is prohibitively difficult. As a result, we are unable to extract 
any insight, from this approach, into the more interesting parts of the phase diagram. 



V. CONCLUSIONS 



In this paper, we report simulation studies of a lattice gas with anisotropic inter-particle interactions, driven to 
a non-equilibrium steady state by an external "electric" field. Focusing first on two-dimensional systems (30 x 30) 
subjected to saturation drive, we find the surprising result that the critical temperature, T c {a, E), is a monotonically 
decreasing function of a = Jn j Jj_ , for fixed J = */ Jji J± (Fig.l). Moreover, it is less than T c (a, 0) of the equilibrium 
system for a > 1.7! These findings remain essentially unchanged for systems up to 90 x 90, leading us to believe that 
T c (a,E) for the infinite system is not likely to be drastically different from those reported here. For drives of the 
order of J, the transitions from the disordered to an ordered phase are continuous. The transition temperatures are 
generally greater than the equilibrium T c 's, at least within the range of a G [1/3,3] explored. The phase diagram is 
more complex (Fig. 2), since two ordered states (a strip aligned along or transverse to the drive) appear, depending on 
a and E. The transition between H- and V-states is, as expected, discontinuous. For reasons yet to be fathomed, the 
line of continuous transitions "dips" at the bi-critical point, where it meets the H-V phase boundary. Finally, we have 
also investigated rectangular systems with aspect ratios which depend on a in such a way that, near criticality, these 
behave isotropically with a simple rescaling of the axes. Though some differences are observed, the phase diagrams 
are qualitatively unchanged. 

On the theoretical front, we have estimated the transition temperatures for large drives, using the simplest method 
available: first order in a "high temperature" expansion. Unfortunately, only the most qualitative feature is repro- 
duced, namely, T c (ao) / 'T C {1 / a®) is greater than unity for ao < 1. Since this method is insensitive to E — O(J), the 
more interesting aspects of the phase diagram are inaccessible. We hope that dynamic mean field methods ]2l}| , which 
can probe small drives, will be successful in producing the behavior in this region of T-a-E space. 

Clearly, this study should be regarded as an initial exploration, providing us with rough guidelines on which parts of 
the phase diagram to probe deeper. Many interesting questions await more detailed simulations as well as theoretical 
analyses. We end with a sample list here. The critical behavior of an equilibrium system is independent of a, i.e., 
whether the system orders into a V- or an H-state. However, for driven systems, renormalized field theory predicts 
drastically different behavior for these two types of ordering j23|. (a) In particular, we expect that transitions into 
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V-states belong to the universality class of the standard driven system, with 5 being the upper critical dimension, etc. 
Thus, checking the validity of this prediction by simulations would be desirable, (b) On the other hand, there are 
no infrared stable fixed points for a theory with longitudinal phase separation, so that this approach fails to describe 
any of the second order transitions into H-states. This kind of ordering is observed in several other non-equilibrium 
systems |^,^6| . An interesting question for simulations studies is whether any of these systems share the same critical 
behavior, thereby establishing new classes of driven diffusive systems, (c) Similar issues arise for the "bi-critical point" . 
For an equilibrium system, due to the underlying symmetry, there can be no new behavior. However, for the driven 
systems, we expect the two critical lines to fall into different universality classes, so that there must be non-trivial 
cross-over phenomena associated with this point. A similar situation arises in the driven bi-layer lattice gas ||. (d) 
In the context of the above issues, finite size scaling [^7j stands out as the common and central analytical tool, (e) 
For the line of first order transitions at lower T, we could improve the accuracy in several ways. This boundary 
was determined by starting runs with a conjectured saddle point configuration ("X"). Checking with other initial 
configurations (e.g., "L") should increase confidence in our findings. Running with more sophisticated algorithms, so 
as to circumvent the long metastable life-times, would be very desirable, (f) Finally, the interfaces in the ordered 
H-states are necessarily inequivalent, since the drive tends to deposit particles on one and remove particles from the 
other. Interfacial correlations in the standard driven systems (ordered in V-states) display remarkable properties p8| . 
A natural question is: what is the nature of the correlation here? In addition, long wavelength interface instabilities 
are quite possible [ p9|]3(]| ], so that it will be extremely interesting to investigate systems with much larger L. From 
this study, it is clear that we are still far from having an intuitive picture of why driven systems order in the way they 
do. Hopefully, with the completion of some of these suggested investigations, we will be closer to the understanding 
of these deceptively simple non-equilibrium systems. 
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